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We present results on the equation of state in QCD with two light quark flavors and a heavier 
strange quark. Calculations with improved staggered fermions have been performed on lattices 
with temporal extent N T = 4 and 6 on a line of constant physics with almost physical quark mass 
values; the pion mass is about 220 MeV, and the strange quark mass is adjusted to its physical 
value. High statistics results on large lattices are obtained for bulk thermodynamic observables, 
i.e. pressure, energy and entropy density, at vanishing quark chemical potential for a wide range 
of temperatures, 140 MeV < T < 800 MeV. We present a detailed discussion of finite cut-off 
effects which become particularly significant for temperatures larger than about twice the transition 
temperature. At these high temperatures we also performed calculations of the trace anomaly on 
lattices with temporal extent N T = 8. Furthermore, we have performed an extensive analysis of zero 
temperature observables including the light and strange quark condensates and the static quark 
potential at zero temperature. These are used to set the temperature scale for thermodynamic 
observables and to calculate renormalized observables that are sensitive to deconfinement and chiral 
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symmetry restoration and become order parameters in the infinite and zero quark mass limits, 
' respectively. 
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I. INTRODUCTION 

Reaching a detailed understanding of bulk thermodynamics of QCD, e.g. the temperature dependence of pressure 
and energy density as well as the equation of state, p(e) vs. e, is one of the central goals of non-perturbative studies of 
■ QCD on the lattice. The equation of state clearly is of central importance for the understanding of thermal properties 
of any thermodynamic system. It provides direct insight into the relevant degrees of freedom and their correlation in 
different phases of strongly interacting matter. We have some understanding of the equation of state in limiting cases 
of high and low temperatures from perturbation theory [1-5] and hadron gas phenomenology [6], respectively. In the 
transition region from the low temperature hadronic regime to the high temperature quark gluon plasma, however, one 
has to rely on a genuine non-perturbative approach, lattice regularized QCD, to study the non-perturbative properties 
of strongly interacting matter. 

Lattice studies of bulk thermodynamics are particularly demanding as the most interesting observables, pressure 
and energy density, are given in terms of differences of dimension 4 operators. These differences are particularly 
difficult to evaluate because both terms being subtracted contain the pressure or energy density of the vacuum, an 
unphysical quantity that is approximately l/(aT) 4 larger than the sought-after difference. Numerical signals thus 
rapidly decrease with the fourth power of the lattice spacing, a, when one tries to approach the continuum limit at 
fixed temperature (T). For this reason improved actions, which allow one to perform calculations on rather coarse 
lattices with relatively small lattice discretization errors, are quite useful in thermodynamic calculations. Indeed, 
the early calculations of bulk thermodynamics with standard staggered [7] and Wilson [8] fermion discretization 
schemes have shown that at high temperature bulk thermodynamic observables are particularly sensitive to lattice 
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discretization errors. This closely follows observations made in studies of the thermodynamics of SU(3) gauge theories 
[9] . In order to minimize discretization errors at high temperature, improved staggered fermion actions - the p4-action 
[10] and the asqtad action [11] - have been used to study the QCD equation of state. Recent studies, performed with 
the asqtad action with almost physical quark mass values on lattices with two different values of the lattice cut-off 
[11], indeed show much smaller discretization errors than similar studies performed with the 1-link, stout smeared 
staggered fermion action [12]. Another source for cut-off errors arises, however, from the explicit breaking of flavor 
symmetry in the staggered fermion formulation. While this is not of much concern in the chirally symmetric high 
temperature phase of QCD, it leads to cut-off dependent modifications of the hadron spectrum and thus may influence 
the calculation of thermodynamic observables in the low temperature hadronic phase of QCD. Techniques to reduce 
flavor symmetry breaking through the introduction of so-called 'fat links' are thus generally exploited in numerical 
calculations with staggered fermions [14]. 

In this paper we report on a calculation of bulk thermodynamics in QCD with almost physical light quark masses 
and a physical value of the strange quark mass. Our calculations have been performed with a tree level Symanzik- 
improved gauge action and an improved staggered fermion action, the p4- action with 3-link smearing (p4fat3), which 
removes 0(a 2 ) cut-off effects at tree-level and also leads to small cut-off effects in 0(g 2 ) perturbation theory [13]. 
At each temperature, we perform simulations with two degenerate light quark masses and a heavier strange quark 
mass for two different values of the lattice cut-off, corresponding to lattices with temporal extent iV r = 4 and 6. 
In these calculations we explore a wide range of temperatures varying from about 140 MeV to about 800 MeV. 
This corresponds to the temperature interval relevant for current experimental studies of dense matter in heavy ion 
collisions at RHIC as well as the forthcoming experiments at the LHC. Bare quark masses have been adjusted to 
keep physical masses approximately constant when the lattice-cut off is varied. At high temperatures, T>350 MeV, 
we also performed calculations on lattices with temporal extent N T — 8 to get control over cut-off effects in the high 
temperature limit. 

We will start in the next section by reviewing basic thermodynamic relations in the continuum valid for thermody- 
namic calculations on such lines of constant physics (LCP). In section III we outline details of our calculational set-up 
with improved staggered fermions. In section IV we present our zero temperature calculations needed to define the 
line of constant physics and the temperature scale deduced from properties of the static quark-antiquark potential. 
Section V is devoted to the presentation of our basic result, the difference between energy density and three times the 
pressure from which we obtain all other thermodynamic observables, e.g. the pressure, energy and entropy densities 
as well as the velocity of sound. Section VI is devoted to a discussion of the temperature dependence of Polyakov loop 
expectation values and chiral condensates which provides a comparison between the deconfining and chiral symmetry 
restoring features of the QCD transition. We finally present a discussion of our results and a comparison with other 
improved staggered fermion calculations of bulk thermodynamics in Section VII. 



II. THERMODYNAMICS ON LINES OF CONSTANT PHYSICS 



To start our discussion of QCD thermodynamics on the lattice we recall some basic thermodynamic relations in the 
continuum. For large, homogeneous media the basic bulk thermodynamic observables we will consider here can be 
obtained from the grand canonical partition function with vanishing quark chemical potentials, Z(T, V). We introduce 
the grand canonical potential, fi(T, V), normalized such that it vanishes at vanishing temperature, 

fi(T, V) = T \nZ{T, V) - Q , (1) 

with £1 = jpiT lnZ(T, V). With this we obtain the thermal part of the pressure (p) and energy density (e) 

,.^.,V,,.4«^, (2) 

which by construction both vanish at vanishing temperature. Using these relations one can express the difference 
between e and 3p, i.e the thermal contribution to the trace of the energy-momentum tensor 6 MM (T), in terms of a 
derivative of the pressure with respect to temperature, i.e. 

(a) 



In fact, it is W (T) which is the basic thermodynamic quantity conveniently calculated on the lattice. All other 
bulk thermodynamic observables, e.g. p/T 4 , e/T 4 as well as the entropy density, s/T 3 = (e + p)/T 4 , can be deduced 
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from this using the above thermodynamic relations. In particular, we obtain the pressure from Qw(T) through 
integration of Eq. 3, 

^-^ = /V^e-(T . (4) 

Usually, the temperature for the lower integration limit, To, is chosen to be a temperature sufficiently deep in the 
hadronic phase of QCD where the pressure p(Tq), receives contributions only from massive hadronic states and is 
already exponentially small. We will discuss this in more detail in Section V. Eq. 4 then directly gives the pressure 
at temperature T. Using p/T 4 determined from Eq. 4 and combining it with Eq. 3, we obtain e/T 4 as well as s/T 3 . 
This makes it evident, that there is indeed only one independent bulk thermodynamic observable calculated in the 
thermodynamic (large volume) limit. All other observables are derived through standard thermodynamic relations so 
that thermodynamic consistency of all bulk thermodynamic observables is insured by construction! 

We stress that the normalization introduced here for the grand canonical potential, Eq. 1, forces the pressure and 
energy density to vanish at T = 0. As a consequence of this normalization, any non-perturbative structure of the 
QCD vacuum, e.g. quark and gluon condensates, that contribute to the trace anomaly 0^(0), and would lead to 
a non- vanishing vacuum pressure and/or energy density, eventually will show up as non-perturbative contributions 
to the high temperature part of these thermodynamic observables. This is similar to the normalization used, e.g. in 
the bag model and the hadron resonance gas, but differs from the normalization used e.g. in rcsummcd pcrturbative 
calculation at high temperature [15, 16] or phenomenological (quasi-particle) models for the high temperature phase of 
QCD [17]. This should be kept in mind when comparing results for the EoS with perturbative and model calculations. 
We also note that ambiguities in normalizing pressure and energy density at zero temperature drop out in a calculation 
of the entropy density which thus is the preferred observable for such comparisons. 



III. LATTICE FORMULATION 



In a lattice calculation, temperature and volume are given in terms of the temporal (N T ) and spatial (N a ) lattice 
extent as well as the lattice spacing, a, which is controlled through the lattice gauge coupling (3 = 6/g 2 , 

T=^y V = (N M P)f- (5) 

As all observables that are calculated on the lattice, are functions of the coupling, (3, we may rewrite Eq. 3 in terms 
of a derivative taken with respect to (3 rather than T . Furthermore we adopt the normalization of the pressure as 
introduced in Eq. 1. This insures a proper renormalization of thermodynamic quantities and, as a consequence, forces 
the pressure to vanish in the vacuum, i.e. at T = 0. 

Let us write the QCD partition function on a lattice of size N^N r as 

Z LCP ((3,N a ,N T )= yn dC/ ^e- S(C/) , (6) 

where U x .^ € SU(3) denotes the gauge link variables and S(U) = (3Sg(U) — Sf(U, 0) is the Euclidean action, which is 
composed out of a purely gluonic contribution, Sg(U), and the fcrmionic part, Sf(U, (3), which arises after integration 
over the fermion fields. We will specify this action in more detail in the next section but note here that we will use 
tree level improved gauge and fermion actions. Although it would be straightforward to introduce one-loop or tadpole 
improvement factors in the action the setup used here greatly simplifies the analysis of thermodynamic observables 
and in some cases also gives a more direct relation to corresponding observables in the continuum. 

When using only tree level improvement the gluonic action does not depend on the gauge coupling, (3, and the 
fermion action depends on (3 only through the bare light (mi) and strange (m s ) quark masses. The subscript LCP 
in Eq. 6 indicates that we have defined the partition function Zlcp on a line of constant physics (LCP), i.e. when 
approaching the continuum limit by increasing the gauge coupling ((3 — > oo) the bare quark masses (rhi(f3), m s ((3)) 
in the QCD Lagrangian are tuned towards zero such that the vacuum properties of QCD remain unchanged. The 
quark masses thus are not independent parameters but are functions of (3 which are determined through constraints 
imposed on zero temperature observables; e.g. one demands that a set of hadron masses remains unchanged when 
the continuum limit is approached on a LCP. 

We now may rewrite Eq. 3 in terms of observables calculable in lattice calculations at zero and non-zero temperature, 

©""CH r. io^( 1 / d5 \ 1 / dS \ \ 
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Here (...) x , with x — r, denote expectation values evaluated on finite temperature lattices of size N^N T , with 
N T -C iV CT , and zero temperature lattices, i.e. on lattices with large temporal extent, N^N T with N T = N >N,j, 
respectively. Furthermore, Rp denotes the lattice version of the QCD /3-function which arises as a multiplicative 
factor in the definition of 9 W (T) because derivatives with respect to T have been converted to derivatives with 
respect to the lattice spacing a on lattices with fixed temporal extent N T , 

We note that in the weak coupling, large (3 limit, Rp approaches the universal form of the 2-loop /3-function of 
3-flavorQCD, 

i? /3 (/?) = 12& o + 72& 1 //? + 0(/r 2 ) , (9) 

with b Q = 9/16tt 2 and b 1 = 1/4tt 4 . 

We analyze the thermodynamics of QCD with two degenerate light quarks (rhi = rh u = rhd) and a heavier strange 
quark (rh s ) described by the QCD partition function given in Eq. 6. For our studies of bulk thermodynamics we use 
the same discretization scheme which has been used recently by us in the study of the QCD transition temperature 
[22], i.e. we use an 0(a 2 ) tree level improved gauge action constructed from a 4-link plaquette term and a planar 
6-link Wilson loop as well as a staggered fermion action that contains a smeared 1-link term and bent 3-link terms. 
We call this action the p4fat3-action; further details are given in Ref. [10] where the p4fat3 action was first used in 
studies of the QCD equation of state on lattices with temporal extent N T — 4 and larger quark masses. With this 
action, bulk thermodynamic quantities like pressure and energy density are 0(a 2 ) tree level improved; corrections to 
the high temperature ideal gas limit only start at 0(1/ N^:) and are significantly smaller than for the Naik action or the 
standard staggered action which suffers from large 0(1/N 2 ) cut-off effects at high temperature. An analysis of cut-off 
effects in the ideal gas limit and in 0(g 2 ) lattice perturbation theory [13] shows that deviations from perturbative 
results are already only a few percent for lattices with temporal extent as small as N T = 6. 

Following the notation used in Ref. [22] the Euclidean action is given as 

S(U)=(3S G (U)-S F (U,(3) , (10) 

with a gluonic contribution, S G (U), and a fermionic part, Sf(U,(3). The latter can be expressed in terms of the 
staggered fermion matrices, £ ) 7 n ! (m s ), for two light (rhi) and a heavier strange quark (rh s ), 

S F (U,f3) = iTrlni?(m J (/3)) + ^Trln^(m s (/3)) . (11) 

Here we took the fourth root of the staggered fermion determinant to represent the contribution of a single fermion 
flavor to the QCD partition function 1 . 

We also introduce the light and strange quark condensates calculated at finite (x = r) and zero temperature (x = 0) , 
respectively, 

W) q ,x = \j^('ftD- 1 {™ q )) x . Q = l,s, x = 0,r, (12) 
as well as expectation values of the gluonic action density, 

(sg)^j^ k (S g ) x . (13) 

All numerical calculations have been performed using the Rational Hybrid Monte Carlo (RHMC) algorithm [21] with 
parameters that have been optimized [22] to reach acceptance rates of about 70%. Some details on our tuning of 
parameters of the RHMC algorithm have been given in [23]. 

For the discussion of the thermodynamics on a line of constant physics (LCP) it sometimes is convenient to 
parametrize the quark mass dependence of Sp in terms of the light quark mass rhi and the ratio h = rh s /rhi rather 
than rhi and rh s separately. We thus write the /3-dependence of the strange quark mass as, rh s (f}) = rhi((3)h((3) . In 



1 There is a controversy regarding the validity of the rooting approach in numerical calculations with staggered fermions. For further 
details we refer to recent reviews presented at Lattice conference [18-20] and references therein. 
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the evaluation of (e — 3p)/T 4 we then will need to know the derivatives of these parametrizations with respect to (3. 
We define 

With these definitions we may rewrite Eq. 7 as 



e ~ 3p = T^(L)=Rm dp/T 



.4 



with 



e{?(T) 



y4 
F 



y4 



r 4 ~ dT Vr 4 / a/3 
egm gro to 



= [<s G >o - <*G>r] JV* , (16) 
= -i^i? m [2m z «^)j, - W>i,r) + m s ((^V) s ,o - W) s> r)] , (17) 
= -R R h m s [<^>,,o - (H>)s,t\ N$ . (18) 



We will show in the next section that to a good approximation h(/3) stays constant on a LCP. Rh thus vanishes on 
the LCP and consequently the last term in Eq. 15, 0^ M , will not contribute to the thermal part of the trace anomaly, 
Q^^(T). The other two terms stay finite in the continuum limit and correspond to the contribution of the thermal 
parts of gluon and quark condensates to the trace anomaly. We note that the latter contribution vanishes in the 
chiral limit of three flavor QCD (rh;, m s — ► 0). The trace anomaly would then entirely be given by 6 G M (T) and 
the observables entering the calculation of bulk thermodynamic quantities in the chiral limit of QCD would reduce 
to those needed also in a pure ST/ (3) gauge theory [9]. In fact, we also find that for physical values of the quark 
masses the trace anomaly is dominated by the gluonic contribution, 6 G M (T). As will become clear in Section V 
B^T) contributes less than 10% to the total trace anomaly for temperatures large than about twice the transition 
temperature. 

We also note that the prefactor in Eq. 17 will approach unity in the continuum limit as R m attains a universal 
form up to 2-loop level which is similar to that of R^ 1 and is only modified through the anomalous dimension of the 
quark mass renormalization [24] . For the relevant combination of j3- functions that enters the fermionic part of the 
trace anomaly, one has 

-ifrWR m Wf-^ = \+™± . (19) 



IV. STATIC QUARK POTENTIAL AND THE LINE OF CONSTANT PHYSICS 
A. Construction of the line of constant physics 

We will calculate thermodynamic observables on a line of constant physics (LCP) that is defined at T — as 
a line in the space of light and strange bare quark masses parametrized by the gauge coupling f3. Each point on 
this line corresponds to identical physical conditions at different values of the lattice cut-off which is tuned towards 
the continuum limit by increasing the gauge coupling (3. We define the line of constant physics by demanding (i) 
that the ratio of masses for the strange pseudo-scalar and the kaon mass, m gs /mK, stays constant and (ii) that rriss 
expressed in units of the scale parameter ro stays constant. The latter gives the distance at which the slope of the zero 
temperature, static quark potential, Vq q (r), attains a certain value. We also introduce the scale r±, which frequently 
is used on finer lattices to convert lattice results expressed in units of the cut-off to physical scales, 

,2 dV qq( r )\ fj2 dV qq( r ) 



=1.65 , =1.0 . (20) 

/ r=r \ / r=n 

We checked that (i) and (ii) also hold true, if we replace m gs by the mass of the light quark pseudo-scalar meson, 
m n . However, errors on m^ro and (m n /m,K) are generally larger which, in particular at large values of (3 makes the 
parametrization of the LCP less stringent. 
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Leading order chiral perturbation theory suggests that the ratio (toss/to^ ) 2 is proportional to rh 8 /(rhi + m 8 ). One 
thus expects this ratio to stay constant for fixed h = to s /to/. This is, indeed fulfilled in the entire regime of couplings, 
(3, explored in our calculations (see Table I) . The first condition for fixing the LCP parameters thus, in practice, has 
been replaced by choosing h = rhi/m s to be constant. As a consequence we find Rh{(3) = 0, which simplifies the 
calculation of thermodynamic quantities. 

In order to define a line of constant strange quark mass, as a second condition for the LCP we demand that the 
product m gs r stays constant. For our LCP we chose 1.59 as the value for the product. Here one should note 
that nigs determined in our calculations only receives contributions from connected diagrams and does not include 
disconnected loops. In order to compare our value (1.59, see discussion below) to a physical one, we therefore follow 
the argumentation of Ref. [25] and adopt m Ss = yj2m 2 K — m\ — 686 MeV as the physical mass of our strange 
pseudoscalar. Together with the scale r = 0.469(7) fm as determined in Ref. [26] through a comparison of r with 
level splittings of the charmonium system [27], this yields m Ss ro ~ 1.63. Of course, there is some ambiguity in this 
choice as current determinations of r differ by about 10% [26, 28]. This introduces some systematic error in the 
definition of the physical LCP. The main reason for deviation from the physical LCP in the present calculation, 
however, is due to the choice of the light quark masses which are about a factor two too large. 

Fixing the light and strange pseudo-scalar masses in units of ro required some trial runs for several (3 values. We 
then used the leading order chiral perturbation theory ansatz m| s ~ m s (or <~ raj) to choose rh s and to/ = to s /10 
at several values of the gauge coupling and used a renormalization group inspired interpolation to determine quark 
mass values at several other (3 values at which high statistics simulations have been performed. It turned out that 
these values are best fitted by m Ss r = 1.59. We thus use this value rather than the value 1.63 mentioned above, to 
define our LCP. For all other simulations we then used the results of these zero temperature calculations to determine 
the quark mass values that belong to a line of constant physics characterized by: 

LCP : (i) m Ss r = 1.59 , (it) h = rh s /rhi = 10 

In general our calculations are thus performed at parameter values close to the LCP which is defined by the above 
condition. The parameters of all our zero temperature calculations performed to determine the LCP, results for 
meson masses and parameters of the static quark potential are summarized in Table I. As can be seen, at our 
actual simulation points the results for m Ss ro fluctuate around the mean value by a few percent. We also checked 
the sensitivity of the meson masses used to determine the LCP to finite volume effects. At (3 = 3.49 and 3.54 we 
performed calculations on 32 4 lattices in addition to the 16 3 • 32 lattices. As can bee seen from Table I results for m gs 
and uiK agree within statistical errors and volume effects are at most on the level of 2% for the light pseudo-scalar. 

The LCP is furthermore characterized by m^/nriK = 0.435(2) and Toss/to^ = 1.33(1). Using ro ~ 0.469(7) fm to 
convert to physical scales we find that on the LCP the light and strange pseudo-scalar masses are ~ 220(4) MeV, 
to 5s — 669(10) MeV and the kaon mass is given by rriK — 503(6) MeV. 

B. The static quark potential and the scale r 

On the LCP we determine several parameters, e.g. the short distance scale r and the linear slope parameter, the 
string tension <r, that characterize the shape of the static quark potential calculated at T = in a fixed range of 
physical distances. The distance ro, defined in Eq. 20, is used to define the temperature scale for the thermodynamics 
calculations. 

The static quark potential, Vq q (r), has been calculated from smeared Wilson loops as described in [22] for all 
parameter sets listed in Table I. We checked that the the smeared Wilson loops project well onto the ground state 
at all values of the cut-off by verifying the independence of the extracted potential parameters on the number of 
smearing levels used in the analysis. The set of gauge couplings, (3 G [3.15,4.08], used in this analysis covers a large 
interval in which the lattice cut-off changes by a factor 6 from a ~ 0.3 fm down to a ~ 0.05 fm. When analyzing 
the static potential over such a wide range of cut-off values one should make sure that the potential is analyzed in 
approximately the same range of physical distances. The fit interval [(r/a) m in, (r/a) m ax] for fits with a Cornell type 
ansatz for the static potential thus has been adjusted for the different values of gauge couplings such that it covers 
approximately the same range of physical distances, r /2 < r < 2r , or 0.25fm<r<lfm. We confirmed our analysis 
of the static quark potential and the determination of ro also independently by using spline interpolations which are 
not biased by a particular ansatz for the form of the potential. 

The left hand part of Fig. 1 shows the static quark potential for several of our parameter sets. We have renormalized 
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6) 


0.3339 


11) 








3.570 


0.212 


24 3 


32 


0.1347 


53) 


0.4053 


18) 


0.3028 


23) 


4.009(26) 


0.2743(38) 


3.414(47) 


3.630 


0.170 


24 3 


32 


0.1126 


20) 


0.3386 


7) 


0.2537 


8) 


4.651(41) 


0.2352(44) 


3.939(59) 


3.690 


0.150 


24 3 


32 


0.1020 


90) 


0.2960 


20) 


0.2230 


30) 


5.201(48) 


0.2116(36) 


4.320(63) 


3.760 


0.130 


24 2 -32-48 


0.0857 


32) 


0.2530 


16) 


0.1894 


16) 


6.050(61) 


0.1810(29) 


4.984(73) 


3.820 


0.125 


32 3 


32 


0.0830 


40) 


0.2310 


38) 


0.1744 


50) 


6.835(44) 


0.1701(21) 


5.541(106) 


3.920 


0.110 


32 3 


32 


0.0750 


70) 


0.2020 


10) 


0.1550 


20) 


7.814(83) 


0.1423(24) 


6.037(72) 


4.080 


0.081 


32 3 


32 


0.0700 


70) 


0.1567 


36) 


0.1220 


50) 


10.39(23) 


0.1060(35) 


7.710(183) 



TABLE I: Light quark and strange pseudo-scalar meson masses and parameters of the static quark potential calculated on zero 
temperature lattices of size N^N T . The last column gives the renormalization constants times ro needed to renormalize the 
heavy quark potential to the string potential at distance r/ro = 1.5. 

these potentials by matching 2 all potentials at a large distance, r/ro = 1.5, to a common value that is taken to be 
identical to the large distance string potential, Turing (>") = — 7r/12r + or. The result of this matching is shown in the 
lower part of Fig. 1 (left) and the constant shifts needed to obtain these renormalized potentials are listed in Table I. 
The good matching of all the potential data obtained at different values of the cut-off already gives a good idea of the 
smallness of finite cut-off effects in this observable. We note that this matching procedure provides renormalization 
constants for the static quark potential, which we also will use later to renormalize the Polyakov loop expectation 
value. 

To further analyze the shape of the static quark potential we determined the scale parameter r /a as well as the 
square root of the string tension in lattice units, *Joa. These parameters have been obtained from three and four 
parameter fits. As described in [22] the latter fit ansatz has been used to estimate systematic errors in our analysis 
of the scale parameters. 

Results for ro/a and y/aa are given in Table I. We note that the product ro^/a stays constant on the LCP and 
changes by less than 2% in the entire range of couplings (3 in which the lattice cut-off changes by a factor 6. For 
a < 0.15 fm we used a quadratic fit ansatz, (ro^/a) a = r^\fo + c(a/ro) 2 , to fit 10 data points. The asymptotic value 
for rQ^/a coincides within errors with a simple average over all values of {ro\/&) a in this interval. This confirms that 
0(a 2 ) corrections indeed are small for this product. Similarly we determined the scale parameter n frequently used 



2 Further details on the matching of the zero temperature heavy quark potentials, its application to the renormalization of the Polyakov 
loop and finite temperature free energies will be published elsewhere. 
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FIG. 1: The static quark potential in units of the scale ro versus distance r/ro (left) and dimensionless combinations of the 
potential shape parameters ro/Vi and roi/a extracted from fits to these potentials (right). The left hand figure shows potentials 
for several values of (3 taken from our entire simulation interval, /3 £ [3.15 : 4.08]. The lowest curve in this figure combines all 
potentials by matching them to the string potential (solid line) as explained in the text. Curves in the right hand figure show 
quadratic fits and a fit to a constant with a 1% error band. The lattice spacing has been converted to physical units using 
r = 0.469 fm. 



to set the scale in calculations performed on finer lattices. Both fits for r^y/a and ro/r\ yield x 2 /dof ~ 0.7. From 
this analysis we obtain the parameters characterizing the shape of the heavy quark potential at masses in the vicinity 
of the LCP, 



r ^ = 1.1034(40) , 
r /ri = 1.4636(60) . 



(21) 



We note that the result obtained here for ro /r± is in good agreement with the corresponding continuum extrapolated 
value, ro/Vi = 1.474(7)(18), determined with the asqtad action from an analysis of the quark mass dependence of this 
ratio at two different values of the lattice spacing, a ~ 0.12 fm and a ~ 0.09 fm, respectively [36]. We show results 
for ro/n and r^^fa calculated at parameter sets close to the LCP in Fig. 1 (right). 

Despite the good scaling behavior of dimensionless combinations of scale parameters deduced from the static 
potential, one expects, of course, to still find substantial deviations from asymptotic scaling relations that are controlled 
by universal 2-loop /^-functions. For the scale parameter r^j a we parametrize deviations from asymptotic scaling using 
a rational function ansatz, 



ro 



a 



l + e r a 2 (/3)+/ r a 4 (/3) 



a r R 2 (p) (1 + b r a 2 (/3) + c r a 4 (/3) + d r a%/3)) 



where 



Ri{P) = exp 



JL 

126r 



-bl/(2b 2 ) 



(22) 



(23) 



denotes the 2-loop /3-function of QCD for three massless quark flavors and a(/3) = i?2(/3)/i?2(3.4). With this 
parametrization it is straightforward to calculate the /3-function Rp entering all basic thermodynamic observables, 



ro 
a 



dr /a 
d/3 



(24) 



Furthermore, we need a parametrization of the /3-dependence of the bare quark masses to determine the second (3- 
function entering the thermodynamic relations, i.e. R m {j3) defined in Eq. 14. For this purpose we use a parametriza- 
tion of the product of the bare light quark mass, rhi and tq that takes into account the anomalous scaling dimension 
of quark masses [24] , 



miro = a m 



126r 



4/9 



PiP) 



(25) 



9 



bm 


Cm 






fm 


9m. 


-2.149(121) 


1.676(178) 


-0.365(144) 


-2.290(162) 


1.829(425) 


-0.356(335) 


a r 


b r 


Cr 




e r 


fr 


13.250(363) 


-1.201(91) 


0.054(196) 


0.406(109) 


-1.682(103) 


0.823(76) 



TABLE II: Parameters of the fit of the scale parameter ro in lattice units based on the ansatz given in Eq. 22 (lower half) and 
the fit of the renormalization group invariant combination of light quark masses and ro (Eq. 25) on the line of constant physics 
(upper half). The x^/d.o.f for these fits are 1.5 for 19 degrees of freedom (lower half) and 0.84 for 18 degrees of freedom (upper 
half). 
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FIG. 2: The scale parameter fo = ro/a versus (3 = 6/g 2 (left) and its product with the bare light quark mass on the LCP 
(right). The two curves shown in the left hand part of this figure correspond to two different fit ansatze. As explained in the 
text in addition to the renormalization group motivated ansatz given in Eq. 22 the result from a 3-interval fit is shown. The 
curve in the right hand part of the figure shows a fit based on the ansatz given in Eqs. 25 and 26. 



with a m being related to the renormalization group invariant quark mass in units of ro and P(/3) being a sixth order 
rational function that parametrizes deviations from the leading order scaling relation for the bare quark mass, 

p(fj) = 1 + b m a 2 (P) + c m a 4 (/3) + d m a 6 (/3) 
W 1 + e m a 2 (P) + f m a\p) + g m a 6 ((3) ' 1 ' 

This ansatz insures that the parametrization for the two /3-functions as well as the parametrization of their product, 
Rp((3)R m ([3), reproduces the universal 2-loop results given in Eqs. 9 and 19. 

In Fig. 2 we show our results for fo = ro/a and m;fo together with the fits described above. The fit parameters 
defining the quark masses on the LCP have been obtained from x 2_n ts in the interval (3 G [3.1,4.08]. Results for the 
fit parameters are given in Table II. In addition we find a m — 0.0190(9) which turns into a value of 8.0(4) MeV in 
physical units. Fit results for ro/a differ from the actually calculated values given in Table II by less than one percent. 

Like in the pure gauge theory calculations of the equation of state, we also find for QCD with light dynamical quarks 
that, in the parameter range of interest for finite temperature calculations, /3-functions deviate significantly from the 
asymptotic scaling form. In particular, we find a dip in Rp at (3 ~ 3.43. For small values of N T , the interesting 
parameter range thus includes the crossover region from the strong to weak coupling regime. 

We use the interpolating fits for fo and m;fo, to determine the two f3- functions Rp and R m that enter the calculations 
of thermodynamic quantities. As all basic thermodynamic observables are directly proportional to Rp, we should check 
the sensitivity of Rp on the particular interpolation form used. We thus have used a completely different interpolation 
that restricts the renormalization group motivated ansatz to the small coupling regime, [3 > 3.52, and uses purely 
rational functions to piecewise fit 2 intervals at smaller (3. We find that results for Rp are sensitive to the fit ansatz 
only for small /3-values, i.e /3<3.25, where the dependence of fo on (3 becomes weak. As discussed later the uncertainty 
on Rp at small values of the coupling only affects the three smallest temperatures used for the analysis of the equation 
of state on the N T — 4 lattices. 
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FIG. 3: The /3-function on the LCP (Eq. 24) (left) and the product RpRm (right). The horizontal lines show the weak coupling 
behavior given in Eqs. 9 and 19. The two curves result from two different fits of fo as discussed in the text. 

Using the parametrizations of fo and rozfo given in Eq. 22 and Eq. 25 as well as the above discussed piecewise 
interpolation of fo we now can derive the two /3-functions Rp((3) and R m {(3). In Fig. 3 we show Rp as well as the 
combination —RpR rri which enter the calculation of the gluonic and fermionic contributions to (e — 3p)/T 4 . For fo 
as well as for the two /3-functions, we show result obtained with our two different fit ansatze. As can be seen, the 
different fit forms lead to differences in the fit result at the edges of the parameter range analyzed. We take care of this 
in our analysis of the equation of state by averaging over the results obtained with the two different fit ansatze and 
by including the difference of both fit results as a systematic error. We note that the /3-function Rp has a minimum 
at (3 ~ 3.43. This characterizes the transition from strong to weak coupling regions and is similar to what is known 
from f3- functions determined in pure gauge theory [9] as well as in QCD with heavier quark masses [10]. The details 
of this region will differ in different discretization schemes as the QCD /3-functions are universal only up to 2-loop 
order in perturbation theory. In order to understand the origin of cut-off effects in thermodynamic observables it is, 
however, important to have good control over Rp in this non-universal regime as well, as Rp enters the calculation of 
all relevant lattice observables as an overall multiplicative factor 3 . 



The basic lattice observables needed to determine the QCD equation of state with our tree level improved gauge and 
fermion actions are expectation values of the gauge action as well as the light and strange quark chiral condensates 
calculated on the LCP on finite (N T <C N a ) and zero (N T >N a ) temperature lattices. We have performed finite 
temperature calculations on lattices with temporal extent N T = 4, 6 and 8. In all cases the spatial extent of the 
lattices was at least four times larger than the temporal extent (N r ), i.e. most finite temperature calculations 
have been performed on lattices of size 16 3 4 and 24 3 6, respectively. In particular at high temperature, we found it 
important to increase the spatial volume in our calculations on N T — 6 lattices to check for possible finite volume 
effects and also to add a few calculations on N T = 8 lattices to get control over the cut-off dependence seen in the 
trace anomaly. In these cases, calculations on 32 3 6 and 32 3 8 lattices have been performed. For all parameter sets, 
corresponding zero temperature calculations have been performed on lattices of size 16 3 32 and 24 3 32. In a few cases 
we used lattices of size 24 2 -32-48 as well as 32 4 . The length of individual calculations on the finite temperature 
lattices varied between 6500 and 35000 trajectories on the N T = 4 lattices and 5000 to 17600 iterations on the N T = 6 
lattices, where Metropolis updates were done after hybrid Monte Carlo evolutions of trajctory length tmd = 0.5. 
At all values of the gauge couplings the length of runs on zero temperature lattices has been adjusted such that the 



We note that in order to insure thermodynamic consistency the /3-function used in the definition of thermodynamic quantities has to 
be determined from the cut-off dependence of the observable used to set the temperature scale, i.e. ro/a in our study. 



V. 



BULK THERMODYNAMICS 



A. The trace anomaly: (e - 3p)/T 4 
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statistical errors of basic observables, e.g. action expectation values, are of similar magnitude as in the T > runs. 
This typically required 2500 to 6000 trajectories. With this amount of statistics, we achieved statistical errors on the 
basic thermodynamic observable, (e — 3p)/T 4 , of below 20% at all temperatures. In fact, they are below 10% in the 
temperature interval T E [180MeV, 700MeV] and are less than 5% for T E [195MeV, 300MeV]. 



0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 




100 200 300 400 500 600 700 800 



FIG. 4: The trace anomaly MAI (T) = e — 3p in units of T 4 versus temperature obtained from calculations on lattices with 
temporal extent 7V T = 4, 6, and 8. The temperature scale, Tro (upper x-axis) has been obtained using the parametrization 
given in Eq. 22, and T [MeV] (lower x-axis), has been extracted from this using ro = 0.469 fm. 

The basic zero and finite temperature observables needed to calculate the trace anomaly in units of the fourth 
power of the temperature, <d^{T)/T 4 = (e - 3p)/T 4 , from Eq. 7 are summarized in Tables III, IV, V and VI. To 
extract MAt (T)/T 4 one furthermore needs to know the derivatives of bare couplings and quark masses, Rp and R m . 
Their calculation from zero temperature observables has been discussed in the previous section. With this input, we 
obtain the result for 9 MA< (T)/r 4 , shown in Fig. 4 for the entire range of temperatures explored by us. Here, and in 
all subsequent figures, the temperature scale has been determined from our results for r^/a, which characterizes the 
slope of the static quark potential and has been extracted from the zero temperature potential as discussed in the 
previous section. On lattices with temporal extent N T we then have Tro = ro/N T . Whenever we show in the following 
temperatures in units of MeV we use r = 0.469 fm [26] to convert Tr to a MeV-scale. We will, however, show in all 
figures both scales which should allow us to compare the results presented here unambiguously with any other lattice 
calculation performed within a different regularization scheme. 

In QCD with light (u, d)-quarks and a heavier strange quark the trace anomaly receives, in addition to the gluonic 
contribution to the trace over the energy-momentum tensor, also contributions from the light and heavy quark chiral 
condensates (Eqs. 16, 17). In the chiral limit only the former contributes and all fermionic contributions enter 
indirectly through modifications of the gauge field background. It thus is interesting to check the relative importance 
of direct contributions from the chiral condensates to (e — 3p)/T 4 . In Fig. 5 we show the fermion contribution Qp^/T 4 
to the total trace anomaly shown in Fig. 4. The right hand part of this figure shows the relative magnitude of the 
light and strange quark contributions. As can be seen they are of similar size close to the transition temperature. 
With increasing temperature, however, the importance of the light quark contribution rapidly drops and becomes 
similar to the ratio of light to strange quark masses at about twice the transition temperature. As can be seen in 
Fig. 5(left) the total fermionic contribution shows a significant cut-off dependence. This partly arises from the large 
change of the product of /3-functions, RpR m that still deviates a lot from the asymptotic weak coupling value in the 
range of couplings relevant for the N T — 4 and 6 calculations, respectively (see Fig. 3(right)). The influence of this 
cut-off dependence on the calculation of the total trace anomaly, however, is strongly reduced as the contribution of 
Qp^/T 4 only amounts to about 20% in the transition region and already drops below 10% at about 1.5T C . 

As all other thermodynamic observables will eventually be deduced from (e — 3p)/T 4 using standard thermodynamic 
relations, we should analyze its structure carefully. Bulk thermodynamics of QCD in different temperature intervals 
addresses quite different physics. This includes (i) the low temperature regime, which in the vicinity of the transition 
temperature often is compared with the physics of a resonance gas and which at lower temperatures is sensitive to 
properties of the hadron spectrum controlled by chiral symmetry breaking; (ii) the genuine non-perturbative physics 
in the transition region and at temperatures above but close to the crossover region which is probed experimentally at 
RHIC and presumably is a still strongly interacting medium with a complicated quasi-particle structure; and (iii) the 
high temperature regime, which eventually becomes accessible to resummed perturbativc calculations. In numerical 



12 



p 


lOOm; 


N* • 




# traj. 


(sg)o 




(tptp)s : o 


3.150 


1.100 
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32 
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32 
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16 3 


32 
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4.31880(34) 
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0.23045(13) 


3.382 
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16 3 


32 
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4.23499(26) 


0.11515(14) 


0.19922(11) 


3.410 


0.412 


16 3 


32 


5824 


4.15990(43) 


0.09013(27) 


0.16256(20) 
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0.390 
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32 
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4.13616(20) 


0.08303(17) 


0.15304(12) 
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0.370 
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0.14364(11) 
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0.344 
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32 
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4.07770(23) 


0.06650(10) 


0.130718(86) 


3.455 


0.329 
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32 
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0.12314(18) 
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32 
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32 3 


32 
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3.396477(51) 
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0.019635(15) 
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0.081 


32 3 
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3.234961(31) 


0.001546(43) 


0.012779(16) 



TABLE III: Expectation values of the pure gauge action density, light and strange quark chiral condensates calculated on zero 
temperature lattices of size N*N T . Also given is the number of trajectories generated at each value of the gauge coupling (3 
with light quarks of mass rhi and bare strange quark mass m s — 10m;. 

calculations on a lattice these, three regimes also deserve a separate discussion as discretization effects influence 
lattice calculations in these regimes quite differently. Before proceeding to a calculation of other bulk thermodynamic 
observables we therefore will discuss in the following three subsections properties of (e — 3p)/T 4 in three temperature 
intervals: (i) T<200 MeV or T<T C , (ii) 200 MeV<T<300 MeV or 1.0<T/T C <1.5 and (hi) T>300 McV or T>1.5T C . 

1. Trace anomaly at low temperatures 

In Fig. 6 we show the low temperature part of (e — 3p) /T 4 obtained from our calculations with the p4fat3 action on 
lattices with temporal extent N T — 4 and 6 and spatial size N a /N T = 4. We compare these results with calculations 
performed with the asqtad action [11] for N T = 6. These latter calculations have been performed on lattices with 
smaller spatial extent, N a /N T = 2, and results are based on lower statistics. These calculations are, however, consistent 
with our findings. We also note that results obtained for two different values of the lattice cut-off, N T — A and 6, are 
compatible with each other. 

In the transition region from high to low temperature it is generally expected that thermodynamic quantities can 
be described quite well by a hadron resonance gas (HRG) [6]; the freeze-out of hadrons in heavy ion experiments 
takes place in this region and observed particle abundances are, in fact, well described by a HRG model [29]. Also a 
comparison of lattice results for the EoS with heavier quarks with a resonance gas model ansatz was quite satisfactory 
[30] but required the use of a suitably adjusted hadron mass spectrum. As we now can perform lattice calculations 
with almost physical quark mass values a more direct comparison using the HRG model with physical quark mass 
values should be appropriate. 

We use a HRG model constructed from all resonances, with masses taken from the particle data book up to a 
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3.74581(27) 


0.0029122(17) 


0.029047(16) 


0.896(49) 


4.4751 


3.690 


0.150 


16 3 


7117 


3.66559(24) 


0.0024312(11) 


0.024276(10) 


0.592(38) 


4.5789 


3.760 


0.130 


16 3 


33378 


3.57727(13) 


0.00199846(36) 


0.0199662(36) 


0.404(22) 


4.6498 


3.820(*) 


0.110 


16 3 


32011 


3.50620(13) 


0.00162776(26) 


0.0162683(26) 


0.273(28) 


4.6830 


3.920 


0.110 


32 3 


6530 


3.395380(89) 


0.00154411(10) 


0.0154337(10) 


0.188(21) 


4.7156 



TABLE IV: Expectation values of the pure gauge action density, light and strange quark chiral condensates calculated on 
lattices with temporal extent N T = 4. The last two columns give the trace anomaly, e — 3p, and the pressure, p, in units of T 4 . 
(*) Note that at j3 = 3.82 simulations on N T = 4 and 6 lattices have been performed at slightly different quark masses. 
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FIG. 5: The fermionic contribution to the trace anomaly (left) and the ratio of the light and strange quark contributions to 
e£7T 4 (right). 
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Here different particle species of mass rrii have degeneracy factors dj and r\i = — 1(+1) for bosons (fermions). A 
comparison of the HRG model with the lattice results is shown as the upper curve in Fig. 6. 

As can be seen in this figure the HRG model captures the qualitative features of the lattice results on (e — 3p) /T 4 
quite well, although the lattice data seem to drop somewhat faster at low temperature. Whether this points towards 
a failure of the HRG model at lower temperatures, or is due to difficulties in correctly resolving the low energy hadron 
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3.335 


0.620 


24 3 


14090 


4.35980(34) 


0.15242(19) 


0.24367(13) 


0.51(25) 


0.0480 


3.351 


0.591 


24 3 


17610 


4.31701(34) 


0.13865(20) 


0.22923(14) 


1.19(25) 


0.0686 


3.382 


0.520 


24 3 


15530 


4.23336(35) 


0.11103(20) 


0.19773(14) 


0.97(19) 


0.1393 


3.410 


0.412 


24 3 


10350 


4.15710(36) 


0.08251(31) 


0.15947(19) 


1.58(24) 


0.2606 


3.420 


0.390 


24 3 


9550 


4.13075(41) 


0.07214(39) 


0.14812(24) 


2.68(19) 


0.3347 


3.430 


0.370 


24 3 


11520 


4.10498(50) 


0.06110(54) 


0.13671(31) 


3.57(23) 


0.4400 


3.445 


0.344 


24 3 


14380 


4.06634(49) 


0.04231(68) 


0.11937(35) 


5.64(28) 


0.6766 


3.455 


0.329 


24 3 


9050 


4.04126(43) 


0.02928(64) 


0.10788(36) 


7.39(32) 


0.8982 


3.460 


0.313 


24 3 


7690 


4.02913(42) 


0.02374(46) 


0.10061(31) 


7.82(34) 


1.0240 


3.470 


0.295 


24 3 


9190 


4.00834(33) 


0.01715(29) 


0.09112(24) 


7.73(26) 


1.2885 


3.490 


0.290 


24 3 


8360 


3.97023(30) 


0.01187(19) 


0.08185(26) 


7.58(27) 


1.7784 


3.510 


0.259 


24 3 


7880 


3.93393(23) 


0.008204(59) 


0.06822(14) 


6.80(20) 


2.2005 


3.540 


0.240 


24 3 


6920 


3.88347(21) 


0.006247(31) 


0.05747(15) 


5.48(23) 


2.7123 


3.570 


0.212 


24 3 


7310 


3.83671(17) 


0.004923(12) 


0.047364(56) 


4.31(19) 


3.0925 


3.630 


0.170 


24 3 


4760 


3.74830(17) 


0.0034263(61) 


0.033892(45) 


2.98(17) 


3.6137 


3.690 


0.150 


24 3 


5190 


3.66697(15) 


0.0027656(24) 


0.027530(19) 


2.09(14) 


3.9362 


3.760 


0.130 


24 3 


8860 


3.57801(12) 


0.0022251(10) 


0.0222031(93) 


1.49(10) 


4.1681 


3.820 


0.125 


32 3 


7870 


3.506568(90) 


0.00203546(42) 


0.0203247(40) 


1.23(11) 


4.3136 


3.920 


0.110 


32 3 


9322 


3.395328(56) 


0.00167642(12) 


0.0167504(12) 


0.973(86) 


4.5057 


4.080 


0.081 


32 3 


6806 


3.234336(54) 


0.00114013(10) 


0.0113976(10) 


0.599(78) 


4.7085 



TABLE V: Expectation values of the pure gauge action density, light and strange quark chiral condensates calculated on lattices 
with temporal extent iV T = 6. The last two columns give the trace anomaly, e — 3p, and the pressure, p, in units of T 4 . 
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(ipip)s,T 
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3.820 
3.920 
4.080 


0.125 
0.110 
0.081 


32 3 
32 3 
32 3 


15100 
27100 
24100 


3.507493(335) 

3.395797(61) 

3.234705(68) 


0.0021449(31) 

0.00174314(39) 

0.00117543(14) 


0.0213654(187) 

0.0174073(31) 

0.0117488(14) 


2.37(23) 
1.72(19) 
0.75(13) 



TABLE VI: Expectation values of the pure gauge action density, light and strange quark chiral condensates calculated on 
lattices with temporal extent N T = 8. The last column gives the trace anomaly, e — 3p, in units of T 4 . 



spectrum in the current calculations on still rather coarse lattices, will require more detailed studies on finer (A^ T = 8) 
lattices in the future. We will return to this question in Section VII. 

We also note that the current lattice calculations are performed with light quark masses that are a factor two larger 
than the physical ones. Reducing the light quark masses to their physical values will shift the lattice data to smaller 
temperatures and will thus improve the comparison with the HRG model. From the known systematics of the quark 
mass dependence of other thermodynamic quantities, e.g. the transition temperature, chiral condensates, or Polyakov 
loop expectation values [22] one can estimate this shift to be less than 5 MeV. Moreover, we note that the scale ro 
used to convert lattice results to physical units has an error of about 2%. This is indicated in Fig. 6 by a horizontal 
error bar for the data. Within this error all data may be shifted coherently. 

The low temperature region of the QCD EoS clearly deserves more detailed study in the future. 

2. Trace anomaly in the strongly non-perturbative regime 

At temperatures just above the transition temperature, (e — 3p)/T 4 shows the largest deviations from the conformal 
limit, e = 3p. The peak in (e— 3p)/T 4 at a temperature T max that is only slightly larger than the transition temperature 
T c constitutes a prominent structure of the trace anomaly which is relatively easy to determine in a lattice calculation. 
It is closely related to the softest point in the QCD equation of state [31], i.e. the minimum of p/e as function of the 
energy density. T max thus plays an important role for the construction of model equations of state that are consistent 
with lattice calculations and may be used in hydrodynamic models for the expansion of dense matter created in heavy 
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FIG. 6: Comparison of the low temperature part of (e — 3p)/T 4 calculated on lattices with temporal extent N T = 4 and 6 with 
a resonance gas model that includes all resonances up to mass 2.5 GeV (dashed curve). The solid curve shows a polynomial fit 
to the N T = 6 data obtained with the p4fat3 action. Data for calculations with the asqtad action are taken from [11]. 



ion collisions. As T max and, in particular (e — 3p)/T^ lax , are easily determined they may also serve as consistency 
checks between different lattice calculations. 

In Fig. 7 we show results for (e — 3p)/T 4 in the intermediate temperature interval 180 MeV < T < 300 MeV. Also 
shown here are results from calculations performed with the asqtad action on lattice with temporal extent N T = 6 
[11]. As can be seen these calculations are in quite good agreement with the results obtained with the p4fat3 action 
on lattice with the same temporal extent but larger spatial volume. Estimates for T max and the peak height on 
lattices with temporal extent N T = 4 and 6 are given in Table VII. Also given there are estimates for the transition 
temperature T c rg obtained previously in a dedicated analysis of the transition temperature on N T = 4 and 6 lattices 
[22]. The values quoted in the Table give the fit results obtained from a joint fit of transition temperatures on both 
lattice sizes for different quark mass values evaluated at the pseudo-scalar mass values corresponding to our LCP. We 
note that the temperature T max is only about 3% larger than the transition temperatures determined from peaks in 
the chiral susceptibility. 

On the coarse N T — 4 lattices the analysis of (e — 3p)/T 4 in the transition region is still quite sensitive to the 
non-perturbative structure of the /3-functions, Rp and RpR m shown in Fig. 3; this region is still close to the strong 
coupling regime below and in the vicinity of the dip in Rp shown in Fig. 3(left). This seems to be the main reason for 
the large differences seen in the peak height for (e — 3p)/T 4 between the N T = 4 and 6 lattices. In the latter case the 
transition and peak region is already in the regime where the lattice /3-functions smoothly approach the continuum 
results. We thus expect that these results are much less affected by this source of lattice artifacts. Nonetheless, a 
better control over the cut-off dependence in this region clearly is needed and does require calculations on a larger 
lattice in order to control the continuum extrapolations of T max ro as well as (e — 3p)/T^ lax . 



3. Trace anomaly at high temperatures 

In Fig. 8 we show results for (e — 3p)/T 4 in the high temperature regime, T>1.5T C . A comparison with data 
obtained with the asqtad action on lattices with temporal extent N T = 6 shows that the results obtained here with 
the p4fat3 action are compatible with the former for T<400 MeV (~ 2T C ). The current analysis performed with the 
p4fat3 action, however, has been extended to much larger temperatures, T ~ 4T C , i.e. into the temperature regime 
accessible to heavy ion experiments at the LHC. 

For temperatures larger than T max the trace anomaly rapidly drops. Eventually, when the high temperature 
perturbative regime is reached, the temperature dependence is expected to be controlled by the logarithmic running 
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FIG. 7: The trace anomaly in the vicinity of the transition temperature calculated on lattices with temporal extent N T = 4 
and 6 on lattices with aspect ratio Ncr/N T — 4. Data for calculations with the asqtad action are taken from [11] which have 
been performed on finite temperature lattices with a smaller physical volume corresponding to an aspect ratio N a /N T — 2. 





T c r 


T c [MeV] 


(TVQ^max 


T max [MeV] 


(£ 


3p) fT max 


4 


0.484(4) 


204(2) 


0.50(1) 


211(4) 




10.8(3) 


6 


0.466(6) 


196(3) 


0.49(1) 


208(4) 




7.8(4) 



TABLE VII: Position of the peak in (e — 3p) /T 4 and its value calculated on lattices with different values of the temporal extent 
iV T on a line of constant physics that corresponds to a pion mass of about 220 MeV. Errors on the peak positions have been 
estimated from cubic fits in the peak region by varying the fit intervals. The second and third columns show the transition 
temperature determined on the LCP used for this study of the EoS. For N T = 4 this had been determined in [22] and for 
N T = 6 in this analysis (see discussion in Section VI). 



of the QCD coupling constant. To leading order in high temperature perturbation theory W (T) for massless quarks 
is given by [1], 

= \b (l + A n/ ) 5 4( T ) + 0(g 5 ) , (28) 

with nf — 3 for massless 3-flavor QCD, which corresponds to the high temperature limit for our (2+l)-flavor QCD 
calculations performed on a LCP with fixed non-zero quark mass values. 

For temperatures larger than about 2.0 T c results for Q^^(T)/T A obviously are sensitive to lattice cut-off effects. 
The results on N T = 6 lattices drop significantly slower with temperature than the N T = 4 results. In order to make 
sure that this effect does not superimpose with possible finite volume effects, we increased in this temperature region 
the spatial lattice size from 24 3 to 32 3 . No statistically significant volume effects have been observed for MM (T)/T 4 , 
although we observe a sensitivity of the zero temperature light and strange quark chiral condensates on the volume; 
as the condensates contribute less than 10% to the trace anomaly at these high temperature values (see Fig. 5) 
modifications of the condensates by a few percent contribute insignificantly to finite volume effects in MM (T)/T 4 . 
Moreover, as the entire fermionic contribution, 0^ I (T)/T 4 , to the total trace anomaly is small for T>400 MeV, it is 
obvious that the contribution of the fermion condensates is not the source for the cut-off effects at high temperature. 
The cut-off dependence seen in Fig. 8 arises from the gluonic sector of 0' i ' i (T)/T 4 , which of course also receives 
contributions from virtual quark loops. 

In the high temperature region we also added calculations on lattices with temporal extent N T = 8 at 3 different 
values of the temperature. Results from these calculations are summarized in Table VI and are also shown in Fig. 8. 
As can be seen in this figure results obtained for the trace anomaly on the N T — 8 lattice are in good agreement with 
the Af T = 6 results suggesting that remaining cut-off effects in this temperature range are small for N T > 6. 

We note that larger values for AtM (T)/T 4 at high temperature also lead to larger values for the pressure, which is 
obtained from an integral over the trace anomaly, and also results in larger values for the energy and entropy densities, 
i.e. these quantities approach the Stefan-Boltzmann limit more rapidly on the iV r = 6 lattices than they did on the 
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FIG. 8: The high temperature part of (e — 3p)/T 4 calculated on lattices with temporal extent N T = 4, 6 and 8. The curves 
show fits to the N T — 4 and 6 data with the ansatz given in Eq. 29. 

N T = 4 lattice. It thus is important to get good control over cut-off effects at high temperatures and obtain further 
confirmation of the results obtained in our N T = 6 lattices, and through further calculations, on the N T — 8 lattices 
at higher temperatures. 

As discussed previously, M/i (T)/T 4 contains a contribution from the vacuum quark and gluon condensates that 
gets suppressed by a factor T 4 at high temperature. In the case of a pure gauge theory it has, however, been noted 
that up to temperatures a few times the transition temperature the dominant power-like correction to the perturbative 
high temperature behavior is 0(T~ 2 ) rather than 0(T~ 4 ) [32, 33]. These qualitative features also show up in our 
results for Q^^{T) /T A at temperatures T>1.5T C . In Fig. 8 we show a comparison of the lattice results with such a 
phenomenologically motivated polynomial fit ansatz, 

/e — 3p\ 3 , a b c 

Here we used the parametric form of the leading order perturbative result given in Eq. 28 with a tempera- 
ture independent coupling g 2 to characterize the high temperature behavior of (e — 3p)/T 4 in the fit interval 
T 6 [300 MeV, 800 MeV] . For N T = 4 we only performed a 2-parameter fit as it turned out that the fit does 
not require the contribution from a constant term (g 2 = 0). The fit parameters obtained from fits in the region 
T > 300 MeV are given in Table VIII. We note that the vacuum condensate contribution (~ c/T 4 ) is small com- 
pared to the genuine thermal part. The present analysis, however, docs not yet allow us to disentangle logarithmic 
from power-like (quadratic) corrections. 

Of course, it is tempting to relate the coefficient of the quartic term to a bag constant or zero temperature quark 
and gluon condensate contribution, c = AB. This yields quite a reasonable value, B 1 / 4 — 247(25) MeV. Nonetheless, 
it seems that a more detailed analysis of the scaling behavior in the high temperature region and better control over 
cut-off effects is needed before a proper running of the gauge coupling can be established that would unambiguously 
allow to single out power-like (quadratic) corrections in the high temperature regime which then would also allow one 
to establish a connection to the perturbative regime for the trace anomaly. 
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g 2 


b [GeV 2 ] 


c [GeV 4 ] 


4 
6 


2.3(7) 


0.101(6) 
0.16(6) 


0.024(1) 
0.013(6) 



TABLE VIII: Fit parameters for 2-parameter fits (N T 
T > 300 MeV using the ansatz given in Eq. 29. 



4) and 3-parameter fits (N T = 6) to (e — 3p)/T 4 in the region 



B. Pressure, energy and entropy density 



As indicated in Eq. 4 we obtain the pressure difference, 

P(T) 



A P (T) 



y4 



P(To) 



(30) 



by integrating over the trace anomaly weighted with an additional factor of T^ 1 in the interval [To,T]. We have 
started our integration at T — 100 MeV, or Tr ~ 0.24, by setting the trace anomaly to zero at this temperature. 
As discussed in the previous section, this leaves us with an uncertainty for the value of the pressure at T , which wc 
estimate to be of the order of the pressure in a hadron resonance gas, i.e. [p(Tq)/Tq]hg = 0.265(2). The results 
obtained for A p (T) from our lattice calculations for the pressure at higher temperatures thus yield p/T 4 up to a 
systematic uncertainty on p(T )/Tq. We also note again that the normalization at T does not take care of the overall 
normalization of the pressure at T = 0. 

To calculate A p (T) by integrating the numerical results obtained for Q^^{T) /T 4 from Eq. 4, we have used straight 
line interpolations of our results for O^/T 4 at adjacent values of the temperature. We also used stepwise interpolations 
obtained by fitting quadratic polynomials to the data in small intervals that are matched to fits in the previous interval. 
Results of the latter approach are then used to perform the integration in the various regions analytically. Differences 
between this approach and the straight line interpolations are nowhere larger than 1.5%. We then used the smooth 
polynomial interpolations to determine the pressure and combined this result with that for 9' i ' i (T) to obtain the 
energy density. Both are shown in the left hand part of Fig. 9. The uncertainty arising from the normalization of 
the pressure at To is indicated as a small vertical bar in the upper right part of this figure. We note that at T ~ 4T C 
results for p/T 4 and e/T 4 stay about 10% below the ideal gas value. 

In particular, for applications to heavy ion phenomenology and for the use of the QCD equation of state in 
hydrodynamic modeling of the expansion of matter formed in heavy ion collisions, it is of importance to eliminate 
the temperature in favor of the energy density and thus obtain the pressure as function of energy density. The ratio 
p/e is shown in the right hand part of Fig. 9. As can be seen at low temperature, in the vicinity of the minimum in 
p/e, results are consistent with values extracted for this quantity from a hadron resonance gas model. We also note 
that in the high temperature regime it has been found in [34] that the ratio p/e shows little dependence the baryon 
number density when evaluated on lines of constant entropy per baryon number. 

The density dependence of p/e is related to the square of the velocity of sound 



dp 



dp/e 
'~de~ 



P 
e 



(31) 



In the high temperature limit as well as in the transition region where the derivative d(p/e)/de vanish, c 2 s is directly 
given by p/e. We therefore find that the velocity of sound is close to the ideal gas value, (? s = 1/3, for energy 
densities e>100 GeV/fm 3 and drops by a factor of 4 to a minimal value of about (c 2 ) m i n ~ 0.09 that is reached at 
e>(l — 2) GeV/fm 3 . The dependence of p/e on the energy density can be parametrized in the high temperature region 
with a simple ansatz [34], 



p - = -\c- 



A 



1 + Be fmVGeV 



(32) 



which then also allows a simple calculation of the velocity of sound, using Eq. 31. We find that the above parametriza- 
tion yields a good fit of the N T = 6 data in the interval 1.3 < e 1 / 4 /(GeV/fm 3 ) 1 / 4 < 6 with a x 2 / dof of 1-3- For the 
fit parameters we obtain, C = 0.964(5), A = 1.16(6) and B = 0.26(3). This fit and the resulting velocity of sound are 
also shown in Fig. 9 (right). 

At energy densities below e ~ 1 GeV/fm 3 the lattice calculations indicate a rise of p/e as expected in hadron 
resonance gas models. However, the current resolution and accuracy of lattice calculations in this regime clearly is 
not yet sufficient to allow for a detailed comparison between both. 
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FIG. 9: Energy density and three times the pressure as function of the temperature (left) and the ratio p/e as function of the 
fourth root of the energy density (right) obtained from calculations on lattices with temporal extent JV T = 4 and 6. Temperature 
and energy density scales have been obtained using the parametrization of ro/a given in Eq. 22 and ro = 0.469 fm. The small 
vertical bar in the left hand figure at high temperatures shows the estimate of the systematic uncertainty on these numbers 
that arises from the normalization of the pressure at To = 100 MeV. The dashed curve (HRG) in the right hand figure shows 
the result for p/e in a hadron resonance gas for temperatures T < 190 MeV. 
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FIG. 10: Entropy density as function of the temperature obtained from calculations on lattices with temporal extent N T — 4 
and 6. Temperature and energy density scales have been obtained using the parametrization of ro/a given in Eq. 22 and 
ro = 0.469 fm. The small vertical bar in the left hand figure at high temperatures shows the estimate of the systematic 
uncertainty on these numbers that arises from the normalization of the pressure at To = 100 MeV. 



As pointed out in Section II the non-perturbative vacuum condensates of QCD show up at high temperature as 
power-like corrections to temperature dependence of the trace anomaly and consequently also to pressure and energy 
density. These vacuum condensate contributions drop out in the entropy density which is shown in Fig. 10. It thus 
is an observable most suitable for comparisons with (resummed) perturbative calculations [15]. Like energy density 
and pressure, the entropy also deviates from the ideal gas value by about 10% at T ~ 4T C . 

We note that for T<2T C the results obtained with the asqtad action [11] for the entropy density are in good 
agreement with the results obtained with the p4fat3 action, although at least in the high temperature limit the cut-off 
dependence of both actions is quite different. This suggests that at least up to temperature T ~ 2T C non-perturbative 
contributions dominate the properties of bulk thermodynamic observables like the entropy density. It also gives rise to 
the expectation that additional cut-off effects are small. Nonetheless, the result presented in this section on properties 
of bulk thermodynamic observables clearly need to be confirmed by calculations on lattices with larger temporal 
extent. 
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VI. RENORMALIZED POLYAKOV LOOP AND CHIRAL CONDENSATES 



As part of our analysis of bulk thermodynamic observables we have gathered a lot of information on the static quark 
potential at zero temperature. This has been discussed in Section IV and results obtained for V qq (r) have been used 
there to determine a temperature scale for our thermodynamic calculations. Furthermore, we have obtained a lot of 
information on the chiral condensates at zero temperature that entered our calculation of thermodynamic quantities. 
Together with corresponding results on heavy quark free energies and chiral condensates at finite temperature this 
allows us to analyze the deconfining properties as well as the change of chiral properties of the finite temperature 
transition in terms of observables which are related to exact order parameters for deconfinement and chiral symmetry 
restoration in the infinite quark mass and vanishing quark mass limits of QCD, respectively. 

As discussed in the previous sections, the deconfining aspect of the finite temperature transition, i.e. the sudden 
liberation of partonic degrees of freedom in QCD, is reflected in the rapid change of bulk thermodynamic observables. 
This is also reflected in the rapid change of the static quark free energy which characterizes the response of a thermal 
medium to the addition of static quark sources. The static quark free energy, F q , is related to the Polyakov loop 
expectation value, (L) ~ exp(— F q (T)/T), 



x = l 



It may more rigorously be defined through the asymptotic large distance behavior of static quark-antiquark correlation 
functions [35], 

(L) 2 = km {L S L\). (34) 

The Polyakov loop needs to be renormalized in order to attain a physically meaningful value in the continuum limit. 
To construct the renormalized Polyakov loop from the bare Polyakov loop expectation values, (L), calculated on 
lattices with temporal extent N T at a temperature controlled by the gauge coupling 0, 

L ren (T) = Z^ n {0){L) , (35) 

we can make use of our extensive calculations of the static potential at zero temperature. As outlined in Section 
IV we have extracted renormalization constants, (c(0)a), from the matching of the static potential to the string 
potential. These renormalization constants are given in Table I in terms of the product c{0)tq. With this we obtain 
the renormalization constants for the Polyakov loop as, Z ren ({3) — exp(c(0)a/2). 

Results for the renormalized Polyakov loop are shown in Fig. 11 (left). We note that the cut-off dependence of L ren 
on lattices with temporal extent N T — 4 and 6 is small, which is in agreement with results obtained in studies of L ren 
in pure ST/ (3) gauge theories [35]. A similar renormalization of the Polyakov loop obtained in calculations with the 
1-link, stout smeared staggered fermion action has been used in [37]. The large cut-off dependence of L ren observed 
in this case mainly seems to be due to the choice of observable {Jk) that has been used to set the temperature scale. 
In fact, when using r$ instead of fx to determine the lattice spacing, and thus the temperature, most of the cut-off 
dependence of L ren is removed in the data shown in [37]. 

Another important aspect of the QCD transition is, of course, the change of chiral properties with temperature. 
This is generally reflected in the temperature dependence of the chiral condensate or related susceptibilities. Also 
the chiral condensates need to be renormalized to obtain finite, well defined quantities in the continuum limit. To 
eliminate the quadratic divergences in the linear quark mass dependent correction to the chiral condensates [24] 
we calculate a suitable combination of light and strange quark condensates at finite temperature. We furthermore 
normalize this quantity by the corresponding combination of condensates calculated at zero temperature at the same 
value of the lattice cut-off, i.e. at the same value of the gauge coupling 0, 

This eliminates multiplicative renormalization factors. 

In Fig. ll(right) we show results obtained for A/. S (T) on the LCP for N T — 4 and 6. In this figure, as well as in the 
corresponding figure for L Ten (T) shown on the left hand side, we also give estimates for the pseudo-critical temperature 
extracted from the position of the peak in the disconnected part of the light quark chiral susceptibility. For N T = 4 
this value has been determined previously by us [22] as the quark mass parameters for the LCP used here are close 
to those used in [22] to determine T c on the N T = 4 lattices for rhi/m s = 0.1. For N T = 6 the choice of the strange 



21 



0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 



0.30 0.40 0.50 0.60 0.70 





100 150 200 250 300 350 400 450 



100 



150 



200 



250 



300 



FIG. 11: Renormalized Polyakov loop on lattices with temporal extent N T — 4, 6 and 8 (left) and the normalized difference of 
light and strange quark chiral condensates defined in Eq. 36. The vertical lines show the location of the transition temperature 
determined in [22] on lattices with temporal extent N T = 4 (right line) and in this analysis for N T = 6 (left line). 



quark mass differs slightly from the one used in that earlier study. We therefore performed a new determination of 
the transition temperature for the N T — 6 lattice and the parameters of the LCP used here. From the peak positions 
of the disconnected parts of the light and strange quark susceptibilities we find (3 C {N T = 6) = 3.445(3). Using the 
value for ro/a quoted for this value of the coupling in Table I we find 4 T c r = 0.466(6) or T c = 196(3). 

We note that the region of most rapid change in the subtracted and normalized chiral condensate, A; iS (T), is in 
good agreement with the region where the Polyakov loop expectation value as well as bulk thermodynamic quantities, 
e.g. the energy and entropy densities change most rapidly. 



VII. DISCUSSION AND CONCLUSIONS 



We have presented here a detailed analysis of the QCD equation of state with an almost physical quark mass 
spectrum. The current calculations have been performed with a physical strange quark mass value and two degenerate 
light quark masses that are about a factor two larger than the physical average light quark mass value. In a wide 
temperature range, results have been obtained on large spatial lattices close to the thermodynamic limit for two 
different values of the lattice cut-off, corresponding to lattices of temporal extent N T — 4 and 6. At high temperatures 
additional calculations on lattices with temporal extent N T = 8 have been performed, which allow us to control 
apparent cut-off effects in this temperature range. All finite temperature calculations have been supplemented with 
corresponding zero temperature calculations to perform necessary vacuum subtractions and to accurately set the 
temperature scale. 

At high temperature, T>2T C , bulk thermodynamic observables such as pressure, energy and entropy density deviate 
from the continuum Stefan-Boltzmann values only by about 10% and show little cut-off dependence. This weak cut- 
off dependence could only be achieved through the use of 0(a 2 ) improved gauge and fermion actions. On the other 
hand, a closer look at the trace anomaly, (e — 3p)/T 4 , from which these quantities are derived, clearly unravels cut-off 
effects when comparing results obtained for the N T = 4 and 6 lattices; for temperatures T>2.5T C or equivalcntly 
T>500 MeV results for (e — 3p)/T 4 on the N T — 4 lattices are systematically lower than for N T — 6. Additional 
calculations performed on N T = 8 lattices in this high temperature region are consistent with the results obtained on 
N T = 6 lattices and thus suggest that cut-off effects are small on lattice with temporal extent N T > 6. Of course, this 
should be confirmed through additional calculations on lattices with temporal extent N T — 8 at larger temperatures. 
On these fine lattices it also will be interesting to analyze in more detail the contribution of charm quarks to the 
equation of state [38, 39]. 



4 Our earlier analysis for m; = 0.1m s on a 16 3 6 lattice, performed with a 20% larger strange quark mass, gave T c ro = 0.4768(51) or 
T c = 201(2). 
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Getting better control over the temperature dependence of (e— 3p)/T 4 at high temperature clearly is important when 
one wants to make contact between lattice calculations for e.g. the entropy density and high temperature perturbation 
theory. Although our present high statistics analysis seems to have achieved good control over the cut-off dependence 
of (e — 3p)/T 4 in this high temperature regime, a more extended analysis of the temperature dependence on N T = 8 
lattices is still needed to make firm contact with perturbative or resummed perturbative calculations. 

At low temperatures, T<200 McV, the influence of cut-off effects is less apparent. We observe that at a given 
value of the temperature results for (e — 3p)/T 4 obtained on the N T = 6 lattice are systematically larger than those 
obtained on the N T = 4 lattice. This is, in fact, expected and is consistent with the cut-off dependence observed in 
calculations of the transition temperature on lattices with temporal extent N T = 4 and 6 [22] . Also on lattices with 
temporal extent N T — 8 indications for this to happen have been found in preliminary studies of chiral and quark 
number susceptibilities as well as Polyakov loop expectation values [40]. We thus expect that with increasing N T , i.e. 
closer to the continuum limit, the region where (e — 3p)/T 4 and all other thermodynamic observables will start to rise 
rapidly, continues to shift towards smaller temperatures. A further, although smaller shift of the transition region 
towards smaller values of the temperature will arise from an extrapolation to physical quark masses. Judging from 
the known temperature dependence of the transition temperature [22] and other thermodynamic observables, like the 
Polyakov loop expectation value or quark number susceptibilities [40] , this will amount to a shift of the scale by a few 
MeV. In fact, extrapolations of the transition temperature in quark mass and lattice spacing to the physical point 
have been performed by several groups for staggered as well as Wilson fermions [22, 41, 42]. These extrapolations 
consistently show that the quark mass dependence of the transition temperature is weak. We take the quark mass 
dependence of the transition temperature as indicator for the shift of the transition region one has to expect in future 
calculations on finer lattices with physical values of the quark masses. We also should stress that current estimates 
for the bare lattice parameters that correspond to physical values of the light quark mass, m q ~ 0.04m s , of course, 
are based on studies of the spectrum on lattices with finite cut-off. Eliminating these systematic effects will require 
further calculations on finer lattices. This may also improve the comparison with model calculations of the equation 
of state in the low temperature phase of QCD. 
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